function output_toolbox_ = dynare_toolbox(M_,options_,order,type,setup_)

global oo_ out
 
oo_.dr.nstatic = M_.nstatic; % this is needed if dynare 4 and above is running

if setup_.order==3 && setup_.pruning==1
oo_  = full_block_dr_new(oo_,M_,options_);
                                           
end

AA   = oo_.dr.ghx ;   
BB   = oo_.dr.ghu ;

if setup_.order >1
    CC   = oo_.dr.ghxx;   DD   = oo_.dr.ghuu;  EE   = oo_.dr.ghxu;   DELTA= oo_.dr.ghs2;  
end

if setup_.order==3 && setup_.pruning==1
    CCC  = oo_.dr.ghxxx; CCD    = oo_.dr.ghxxu; DDD    = oo_.dr.ghuuu; EED  = oo_.dr.ghxuu; DELTAA = oo_.dr.ghs2x; DELTAB = oo_.dr.ghs2u;
end

ys                = oo_.dr.ys(oo_.dr.order_var,:);
exosize           = M_.exo_nbr;   
select_state      = M_.nstatic+1 : M_.nstatic + M_.npred + M_.nboth;
M_.nspred         = M_.npred + M_.nboth;
endosize          = length(oo_.dr.ghx);                                 
statesize         = M_.nspred;
sim_first         = zeros(M_.endo_nbr,setup_.lengthsimul);
sim_second        = zeros(M_.endo_nbr,setup_.lengthsimul);
sim_third         = zeros(M_.endo_nbr,setup_.lengthsimul);
sim_first_sigma_2 = zeros(M_.endo_nbr,setup_.lengthsimul);
sim_total         = zeros(M_.endo_nbr,setup_.lengthsimul);
simulstore        = zeros(setup_.numbersimul*endosize,setup_.lengthsimul);
shocks_matrix     = zeros(setup_.numbersimul*exosize,setup_.lengthsimul); 


if setup_.hpfilter==1 || setup_.gen_irf==1
    
    for k_index = 1:exosize
        for s_index   = 1:setup_.numbersimul
            rng(142857*k_index + s_index); %this ensures you always draw the same draws for each shock
            ndraws     = randn(1,setup_.lengthsimul);
          
            indeces = find(abs(ndraws)>2); 

            if size(indeces)>0
               for h_index=1:size(indeces,2)
                   while abs(ndraws(indeces(h_index)))>2
                       ndraws(indeces(h_index))  = randn(1,1);
                   end
               end
            end
          
            epsZ                                             = setup_.stdev_shocks(:,k_index)*ndraws;
            shocks_matrix(k_index + (s_index - 1)*exosize,:) = epsZ;
            
            if (options_.periods~=0)
                shocks_matrix     = oo_.exo_simul';
            end
            
        end
    end

    if setup_.gen_irf==1
            shocksM_half1 = shocks_matrix(1:1:exosize*setup_.numbersimul/2,:); % take first half of shocks_matrix
            shocksM_half2 = shocksM_half1;                                     % take second half and set equal to first half
            for k_index = 1:exosize                                            % change epsZ at period setup_.drop for shocks that you consider for IRF (in each simulation)
                if setup_.stdev_shocks(1,k_index)~=0
                    for s_index = 1: setup_.numbersimul/2
                        shocksM_half2(k_index + (s_index - 1)*exosize,setup_.drop) = shocksM_half2(k_index + (s_index - 1)*exosize,setup_.drop)+ setup_.stdev_shocks(1,k_index);
                    end
                end
            end
            shocks_matrix(exosize*setup_.numbersimul/2 + 1:1:exosize*setup_.numbersimul,:) = shocksM_half2;
            checkM = max((max(abs(shocksM_half1 - shocksM_half2)).^2)); max_effective_stdev = max(setup_.stdev_shocks.*setup_.stdev_shocks);
            if (abs(checkM - max_effective_stdev)>1e-15) || checkM==0; disp('error shocks structure with setup_.gen_irf==1: dynare_toolbox'); open dynare_toolbox; keyboard; end 
    end
end


if setup_.stoch_ss_irf==1
    setup_.numbersimul       = 1;
    shocks_matrix            = zeros(exosize,setup_.lengthsimul);
    simulstore               = zeros(setup_.numbersimul*endosize,setup_.lengthsimul);

    
       shocks_matrix(1, setup_.drop)      =   setup_.stdev_shocks(:,1);   
       shocks_matrix(2, setup_.drop)      =   setup_.stdev_shocks(:,2);   
       
    
    shocks_matrix(1,:) = setup_.scale_shock_1*shocks_matrix(1,:);

    if (options_.periods~=0)
        shocks_matrix     = oo_.exo_simul';
    end
end

if setup_.external_shocks==1
      shocks_matrix = setup_.shocks_matrix_external;
end


for s_index = 1:setup_.numbersimul
    initialValues   = ys(select_state,1); 
    statesInitial   = initialValues - ys(select_state,1);
    Shocks          = shocks_matrix(1 + (s_index-1)*exosize: exosize*s_index, :); 
    
    if setup_.order == 1
        sim_first(:,1)  = AA*statesInitial + BB*Shocks(:,1);
        sim_total(:,1)  = ys + sim_first(:,1);
        for t=2:setup_.lengthsimul
            sim_first(:,t) = AA*sim_first(select_state,t-1)+BB*Shocks(:,t);
            sim_total(:,t) = ys + sim_first(:,t);
        end
    end
    
    if setup_.pruning == 0 && setup_.order>1
        
        sim_first(:,1)  = AA*statesInitial + BB*Shocks(:,1);
        kronShocks      = kron(Shocks(:,1),Shocks(:,1)); kronXhShock = kron(statesInitial,Shocks(:,1)); kronXh = kron(statesInitial,statesInitial);
        sim_second(:,1) = (1/2)*(CC*kronXh+2*EE*kronXhShock+DD*kronShocks+DELTA);  
        
        if setup_.order==2
            sim_total(:,1)  = ys + sim_first(:,1) + sim_second(:,1);            
            for t=2:setup_.lengthsimul
                sim_first(:,t)  = AA*sim_first(select_state,t-1)+BB*Shocks(:,t);
                kronShocks      = kron(Shocks(:,t),Shocks(:,t)); kronXhShock = kron((sim_total(select_state,t-1) - ys(select_state)),Shocks(:,t)); kronXh = kron((sim_total(select_state,t-1) - ys(select_state)),(sim_total(select_state,t-1) - ys(select_state)));
                sim_second(:,t) = AA*sim_second(select_state,t-1)+(1/2)*(CC*kronXh+2*EE*kronXhShock+DD*kronShocks+DELTA);
                sim_total(:,t)  = ys + sim_first(:,t) + sim_second(:,t);
            end
        end
        
        if setup_.order==3
            n                       = M_.nspred + M_.exo_nbr; 
            G0                      = oo_.dr.g_0;   % matrix rows correspond to all endogenous in DR-order
            G1                      = oo_.dr.g_1;   % matrix rows correspond to all endogenous in DR-order
            G2                      = zeros(endosize,n*n);
            G3                      = zeros(endosize,n*n*n);

            i = 1;
            while i < n+1
                eval(['oo_.dr.g_2' num2str(i) '= oo_.dr.g_2(:,(- n + (n + 1.5)*i - (i^2)/2): (i*n - (0.5*(i^2) - (0.5*i))) );'])
                i = i+1;
            end
            
            for i = 1:n
                G2(:,i*(n+1)-n:i*n) = eval(['oo_.dr.g_2' num2str(i)]);  % position the partitions in their respective places
                for k = 2:n
                    for j = k:i
                        G2(:,(j-(k-1))*n+(k-2)*n+(k-1)) = eval(['oo_.dr.g_2' num2str(k-1) '(:,j-(k-2));']);  % position the partitions elements (columns) to set the symmetric terms
                    end
                end
            end
            
            m = 1;
            for i=1:n
                for j=i:n
                    for k=j:n
                        xx = oo_.dr.g_3(:,m);
                        G3(:,(i-1)*n*n+(j-1)*n+k) = xx;
                        G3(:,(i-1)*n*n+(k-1)*n+j) = xx;
                        G3(:,(j-1)*n*n+(k-1)*n+i) = xx;
                        G3(:,(j-1)*n*n+(i-1)*n+k) = xx;
                        G3(:,(k-1)*n*n+(i-1)*n+j) = xx;
                        G3(:,(k-1)*n*n+(j-1)*n+i) = xx;
                        m = m + 1;
                    end
                end
            end
       
        z_dyn                   = zeros(statesize + exosize,setup_.lengthsimul);       % matrix of state variables + exo shocks
        z_dyn(:,1)              = vertcat(statesInitial,Shocks(:,1));
        sim_total(:,1)          = ys + G0 + G1*z_dyn(:,1) + G2*kron(z_dyn(:,1),z_dyn(:,1)) + G3*kron(kron(z_dyn(:,1),z_dyn(:,1)),z_dyn(:,1));
        for t=2:setup_.lengthsimul
            z_dyn(:,t)          = vertcat(sim_total(select_state,t-1) - ys(select_state), Shocks(:,t));
            sim_total(:,t)      = ys + G0 + G1*z_dyn(:,t) + G2*kron(z_dyn(:,t),z_dyn(:,t)) + G3*kron(kron(z_dyn(:,t),z_dyn(:,t)),z_dyn(:,t));     
        end
        end
        
        if setup_.order==4

            create_g_matrix;

            z_dyn                   = zeros(statesize + exosize,setup_.lengthsimul);       % matrix of state variables + exo shocks
            z_dyn(:,1)              = vertcat(zeros(statesize,1),Shocks(:,1));
            z_z                     = kron(z_dyn(:,1),z_dyn(:,1));
            z_z_z                   = kron(kron(z_dyn(:,1),z_dyn(:,1)),z_dyn(:,1));
            z_z_z_z                 = kron(z_dyn(:,1),z_z_z);
            sim_total(:,1)          = ys + oo_.g0_dp + (oo_.g1_dp + oo_.g1ssx_dp + oo_.g1sssx_dp)*z_dyn(:,1) + (oo_.g2_dp + oo_.g2ssxx_dp)*z_z + oo_.g3_dp*z_z_z + oo_.g4_dp*z_z_z_z;


            for t=2:setup_.lengthsimul
                z_dyn(:,t)          = vertcat(sim_total(select_state,t-1) - ys(select_state), Shocks(:,t));
                z_z                 = kron(z_dyn(:,t),z_dyn(:,t));
                z_z_z               = kron(kron(z_dyn(:,t),z_dyn(:,t)),z_dyn(:,t));
                z_z_z_z             = kron(z_dyn(:,t),z_z_z); 
                sim_total(:,t)      = ys + oo_.g0_dp + (oo_.g1_dp + oo_.g1ssx_dp + oo_.g1sssx_dp)*z_dyn(:,t) + (oo_.g2_dp+oo_.g2ssxx_dp)*z_z + oo_.g3_dp*z_z_z + oo_.g4_dp*z_z_z_z;
            end 
        end
    end

    
    if setup_.pruning==1
        
        if setup_.order==2
            
            sim_first(:,1)  = AA*statesInitial + BB*Shocks(:,1);
            kronShocks      = kron(Shocks(:,1),Shocks(:,1)); kronXhShock = kron(statesInitial,Shocks(:,1)); kronXh = kron(statesInitial,statesInitial);
            
            if strcmp(type,'kim_et_al')
                sim_second(:,1) = (1/2)*(CC*kronXh+2*EE*kronXhShock+DD*kronShocks+DELTA);
                sim_total(:,1)  = ys + sim_first(:,1) + sim_second(:,1);      
                for t=2:setup_.lengthsimul
                    sim_first(:,t)  = AA*sim_first(select_state,t-1)+BB*Shocks(:,t);
                    kronShocks      = kron(Shocks(:,t),Shocks(:,t)); kronXhShock = kron(sim_first(select_state,t-1),Shocks(:,t)); kronXh = kron(sim_first(select_state,t-1),sim_first(select_state,t-1));
                    sim_second(:,t) = AA*sim_second(select_state,t-1)+(1/2)*(CC*kronXh+2*EE*kronXhShock+DD*kronShocks+DELTA);
                    sim_total(:,t)  = ys + sim_first(:,t) + sim_second(:,t);
                end
            end
            if strcmp(type,'den_haan_de_wind')
                sim_second(:,1) = (1/2)*(CC*kronXh+2*EE*kronXhShock+DD*kronShocks);
                sim_total(:,1)  = ys + (1/2)*DELTA + sim_first(:,1) + sim_second(:,1);      
                for t=2:setup_.lengthsimul
                    sim_first(:,t)  = AA*sim_first(select_state,t-1)+BB*Shocks(:,t);
                    kronShocks      = kron(Shocks(:,t),Shocks(:,t));kronXhShock=kron(sim_first(select_state,t-1),Shocks(:,t));kronXh=kron(sim_first(select_state,t-1),sim_first(select_state,t-1));
                    sim_second(:,t) = AA*sim_second(select_state,t-1)+ (1/2)*(CC*kronXh+2*EE*kronXhShock+DD*kronShocks);
                    sim_total(:,t)  = ys + (1/2)*DELTA + sim_first(:,t) + sim_second(:,t);
                end
            end
        end
        
        if setup_.order==3
            sim_first(:,1)         =  AA*statesInitial + BB*Shocks(:,1);
            kronShocks             = kron(Shocks(:,1),Shocks(:,1)); kronXhShock = kron(statesInitial,Shocks(:,1)); kronXh = kron(statesInitial,statesInitial);
            sim_second(:,1)        = (1/2)*(CC*kronXh+2*EE*kronXhShock+DD*kronShocks+DELTA);
            sim_first_sigma_2(:,1) = (1/2)*(DELTAA*statesInitial+DELTAB*Shocks(:,1));
            sim_third(:,1)         = (1/6)*(CCC*kron(statesInitial,kronXh)+DDD*kron(Shocks(:,1),kronShocks))...
                                     +(1/2)*(CCD*kron(kronXh,Shocks(:,1))+EED*kron(statesInitial,kronShocks))...
                                     + CC*kron(statesInitial,statesInitial)+EE*kron(statesInitial,Shocks(:,1));
            sim_total(:,1)         = ys + sim_first(:,1) + sim_second(:,1) + sim_third(:,1) + sim_first_sigma_2(:,1); 
            
            if strcmp(type,'andreasen') % up to second order identical to kim_et_al
                for t=2:setup_.lengthsimul
                    sim_first(:,t)         = AA*sim_first(select_state,t-1)+BB*Shocks(:,t);
                    kronShocks             = kron(Shocks(:,t),Shocks(:,t)); kronXhShock = kron(sim_first(select_state,t-1),Shocks(:,t)); kronXh = kron(sim_first(select_state,t-1),sim_first(select_state,t-1));
                    sim_second(:,t)        = AA*sim_second(select_state,t-1) +(1/2)*(CC*kronXh +2*EE*kronXhShock  +DD*kronShocks +DELTA);
                    sim_first_sigma_2(:,t) = AA*sim_first_sigma_2(select_state,t-1)+(1/2)*(DELTAA*sim_first(select_state,t-1)+DELTAB*Shocks(:,t));
                    sim_third(:,t)         = AA*sim_third(select_state,t-1)...
                                             +(1/6)*(CCC*kron(sim_first(select_state,t-1),kronXh)+DDD*kron(Shocks(:,t),kronShocks))...
                                             +(1/2)*(CCD*kron(kronXh,Shocks(:,t))+EED*kron(sim_first(select_state,t-1),kronShocks))...
                                             + CC*kron(sim_second(select_state,t-1),sim_first(select_state,t-1))...
                                             + EE*kron(sim_second(select_state,t-1),Shocks(:,t));
                    sim_total(:,t)         = ys + sim_first(:,t) + sim_second(:,t) + sim_third(:,t) + sim_first_sigma_2(:,t);
                end
            end
            
            if strcmp(type,'fernandez-villaverde_et_al')% up to second order identical to kim_et_al
                for t=2:setup_.lengthsimul
                    sim_first(:,t)         = AA*sim_first(select_state,t-1)+BB*Shocks(:,t);
                    kronShocks             = kron(Shocks(:,t),Shocks(:,t)); kronXhShock=kron(sim_first(select_state,t-1),Shocks(:,t)); kronXh=kron(sim_first(select_state,t-1),sim_first(select_state,t-1));
                    sim_second(:,t)        = AA*sim_second(select_state,t-1)+(1/2)*(CC*kronXh+2*EE*kronXhShock+DD*kronShocks+DELTA);
                    sim_first_sigma_2(:,t) = AA*sim_first_sigma_2(select_state,t-1)+(1/2)*(DELTAA*sim_first(select_state,t-1)+DELTAB*Shocks(:,t));
                    sim_third(:,t)         = AA*sim_third(select_state,t-1)...
                                             +(1/6)*(CCC*kron(sim_first(select_state,t-1),kronXh)+DDD*kron(Shocks(:,t),kronShocks))...
                                             +(1/2)*(CCD*kron(kronXh,Shocks(:,t))+EED*kron(sim_first(select_state,t-1),kronShocks));
                    sim_total(:,t)         = ys + sim_first(:,t) + sim_second(:,t) + sim_third(:,t) + sim_first_sigma_2(:,t);                    
                end
            end
            
            if strcmp(type,'juillard')% up to second order identical to kim_et_al
                for t=2:setup_.lengthsimul
                    sim_first(:,t)         = AA*sim_first(select_state,t-1)+BB*Shocks(:,t);
                    kronShocks             = kron(Shocks(:,t),Shocks(:,t)); kronXhShock=kron(sim_first(select_state,t-1),Shocks(:,t)); kronXh=kron(sim_first(select_state,t-1),sim_first(select_state,t-1));
                    sim_second(:,t)        = AA*sim_second(select_state,t-1)+(1/2)*(CC*kronXh+2*EE*kronXhShock+DD*kronShocks+DELTA);
                    sim_first_sigma_2(:,t) = AA*sim_first_sigma_2(select_state,t-1)+(1/2)*(DELTAA*sim_first(select_state,t-1)+DELTAB*Shocks(:,t));
                    sim_third(:,t)         = AA*sim_third(select_state,t-1)...
                                             +(1/6)*(CCC*kron(sim_first(select_state,t-1),kronXh)+DDD*kron(Shocks(:,t),kronShocks))...
                                             +(1/2)*(CCD*kron(kronXh,Shocks(:,t))+EED*kron(sim_first(select_state,t-1),kronShocks))...
                                             +CC*kron(sim_second(select_state,t-1),sim_first(select_state,t-1));
                   sim_total(:,t)         = ys + sim_first(:,t) + sim_second(:,t) + sim_third(:,t) + sim_first_sigma_2(:,t); 
                end
            end
            
            if strcmp(type,'den_haan_de_wind')
                sim_first(:,1)         = (AA+(1/2)*DELTAA)*statesInitial+(BB+(1/2)*DELTAB)*Shocks(:,1);
                kronShocks             = kron(Shocks(:,1),Shocks(:,1));kronXhShock=kron(statesInitial,Shocks(:,1));kronXh=kron(statesInitial,statesInitial);
                sim_second(:,1)        = (1/2)*DELTAA*statesInitial+(1/2)*(CC*kronXh+2*EE*kronXhShock+DD*kronShocks);
                sim_third(:,1)         = (1/2)*DELTAA*statesInitial+(BB+(1/2)*DELTAB)*Shocks(:,1)...
                                         +(1/2)*(CC*kron(statesInitial,statesInitial)+2*EE*kron(statesInitial,Shocks(:,1))+DD*kronShocks)...
                                         +(1/6)*(CCC*kron(statesInitial,kronXh)+DDD*kron(Shocks(:,1),kronShocks))...
                                         +(1/2)*(CCD*kron(kronXh,Shocks(:,1))+EED*kron(statesInitial,kronShocks));
               sim_total(:,1)          = ys + (1/2)*DELTA + sim_third(:,1);
               for t=2:setup_.lengthsimul
                   sim_first(:,t)      = (AA+(1/2)*DELTAA)*sim_first(select_state,t-1)+(BB+(1/2)*DELTAB)*Shocks(:,t);
                   kronShocks          = kron(Shocks(:,t),Shocks(:,t));kronXhShock=kron(sim_first(select_state,t-1),Shocks(:,t));kronXh=kron(sim_first(select_state,t-1),sim_first(select_state,t-1));
                   sim_second(:,t)     = (AA+(1/2)*DELTAA)*sim_second(select_state,t-1)...
                                         +(1/2)*(CC*kronXh+2*EE*kronXhShock+DD*kronShocks);
                   sim_third(:,t)      = (AA+(1/2)*DELTAA)*sim_third(select_state,t-1)+(BB+(1/2)*DELTAB)*Shocks(:,t)...
                                         +(1/2)*(CC*kron(sim_second(select_state,t-1),sim_second(select_state,t-1))+2*EE*kron(sim_second(select_state,t-1),Shocks(:,t))+DD*kronShocks)...
                                         +(1/6)*(CCC*kron(sim_first(select_state,t-1),kronXh)+DDD*kron(Shocks(:,t),kronShocks))...
                                         +(1/2)*(CCD*kron(kronXh,Shocks(:,t))+EED*kron(sim_first(select_state,t-1),kronShocks));
                   sim_total(:,t)      = ys + (1/2)*DELTA + sim_third(:,t);
               end
            end
            
            if strcmp(type, 'lan_meyer-gohde')
              global abounds_; %#ok<TLEV>
              abounds_.Shocks           = Shocks; %#ok<*STRNU>
              abounds_.lengthsimul      = setup_.lengthsimul;
              pruning_abounds(M_,options_,order,type);
              sim_total                 = abounds_.sim_total; 
            end

        end
        
        if setup_.order==4  

            shocks                 = Shocks;
            sig                    = 1;
            out.eta(out.eta~=0)    = 1;

                        
            ny      = size(out.gv,1);        %Number of output variables
            nx      = size(out.hv,1);        %Number of state variables
            
 
            %Defining matrices
            HHxxtil   = 1/2*reshape(out.hvv,nx,nx^2);
            HHxxxtil  = 1/6*reshape(out.hvvv,nx,nx^3);
            HHxxxxtil = 1/24*reshape(out.hvvvv,nx,nx^4);
            HHssxxtil = 6/24*reshape(out.hssvv,nx,nx^2);
            GGxxtil   = 1/2*reshape(out.gvv,ny,nx^2);   
            GGxxxtil  = 1/6*reshape(out.gvvv,ny,nx^3);   
            GGxxxxtil = 1/24*reshape(out.gvvvv,ny,nx^4);   
            GGssxxtil = 6/24*reshape(out.gssvv,ny,nx^2);   
            
            %% FIND DEVIATIONS FROM NON-STOCHASTIC STEADY STATE
            xf  = zeros(nx,1);
            xs  = ((-  out.hv + eye(size(out.hv)))^-1)*1/2*sig^2*out.hss;  
            xrd = ((-  out.hv + eye(size(out.hv)))^-1)*1/6*sig^3*out.hsss;
            xth = ((-  out.hv + eye(size(out.hv)))^-1)*(HHxxtil*( kron(xs,xs)) + 3/6*out.hssv*sig^2*xs + 1/24*out.hssss*sig^4);
            
  
            Y_sim     = zeros(ny,setup_.lengthsimul);
            X_sim     = zeros(nx,setup_.lengthsimul);
            
            Y_sim_1st   = 0*Y_sim;
            Y_sim_2nd   = 0*Y_sim;
            Y_sim_3rd   = 0*Y_sim;
            Y_sim_4th   = 0*Y_sim;
  
            X_sim_1st   = 0*X_sim;
            X_sim_2nd   = 0*X_sim;
            X_sim_3rd   = 0*X_sim;
            X_sim_4th   = 0*X_sim;

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            %SET INITIAL STATES TO STOCHASTIC STEADY STATE
            xf_xf       = kron(xf,xf);
            xf_xs       = kron(xf,xs);
            xf_xf_xf    = kron(xf,xf_xf);
            xf_xf_xf_xf = kron(xf,xf_xf_xf);
            xf_xrd      = kron(xf,xrd);
            xf_xf_xs    = kron(xf_xf,xs);
            xs_xs       = kron(xs,xs);
            xf_p        = out.hv*xf  + sig*out.eta*shocks(:,1);
            xs_p        = out.hv*xs  + HHxxtil*xf_xf   + 1/2*sig^2*out.hss;
            xrd_p       = out.hv*xrd + 2*HHxxtil*xf_xs + HHxxxtil*xf_xf_xf + 3/6*out.hssv*sig^2*xf + 1/6*sig^3*out.hsss;
            xth_p       = out.hv*xth + 2*HHxxtil*xf_xrd + HHxxtil*xs_xs + 3*HHxxxtil*xf_xf_xs + HHxxxxtil*xf_xf_xf_xf + 3/6*out.hssv*sig^2*xs + HHssxxtil*sig^2*xf_xf + 4/24*out.hsssv*sig^3*xf + 1/24*out.hssss*sig^4;
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            %needed at order 2
            xf_xf_p    = kron(xf_p,xf_p);

            %needed at order 3
            xf_xs_p    = kron(xf_p,xs_p);
            xf_xf_xf_p = kron(xf_p,xf_xf_p);

            %needed at order 4
            xf_xf_xf_xf_p = kron(xf_p,xf_xf_xf_p);
            xf_xrd_p      = kron(xf_p,xrd_p);
            xf_xf_xs_p    = kron(xf_xf_p,xs_p);
            xs_xs_p       = kron(xs_p,xs_p);

   
            %CONTROL VARIABLES
            Y_sim_1st(:,1) = out.gv*xf_p;
            Y_sim_2nd(:,1) = out.gv*xs_p  + GGxxtil*xf_xf_p + 1/2*sig^2*out.gss;
            Y_sim_3rd(:,1) = out.gv*xrd_p + GGxxtil*2*xf_xs_p + GGxxxtil*xf_xf_xf_p + 3/6*out.gssv*sig^2*xf_p + 1/6*sig^3*out.gsss;       
            Y_sim_4th(:,1) = out.gv*xth_p + GGxxtil*(2*xf_xrd_p+xs_xs_p) + GGxxxtil*(2*xf_xf_xs_p) +...
                             GGxxxxtil*xf_xf_xf_xf_p + 3/6*out.gssv*sig^2*xs_p + GGssxxtil*sig^2*xf_xf_p + 4/24*out.gsssv*sig^3*xf_p + 1/24*sig^4*out.gssss;
   
            %UPDATE
            xf       = xf_p;
            xs       = xs_p;        
            xrd      = xrd_p;
            xth      = xth_p;

            %order 2
            xf_xf    = xf_xf_p;

            %order 3
            xf_xs    = xf_xs_p;
            xf_xf_xf = xf_xf_xf_p;

            %order 4
            xf_xf_xf_xf = xf_xf_xf_xf_p;
            xf_xrd      = xf_xrd_p;
            xf_xf_xs    = xf_xf_xs_p;
            xs_xs       = xs_xs_p;

            %RUN LOOP
            for t=2:setup_.lengthsimul

                xf_p              = out.hv*xf  + sig*out.eta*shocks(:,t);
                xs_p              = out.hv*xs  + HHxxtil*xf_xf   + 1/2*sig^2*out.hss;
                xrd_p             = out.hv*xrd + 2*HHxxtil*xf_xs + HHxxxtil*xf_xf_xf + 3/6*out.hssv*sig^2*xf + 1/6*sig^3*out.hsss;
                xth_p             = out.hv*xth + 2*HHxxtil*xf_xrd + HHxxtil*xs_xs + 3*HHxxxtil*xf_xf_xs + HHxxxxtil*xf_xf_xf_xf + 3/6*out.hssv*sig^2*xs + HHssxxtil*sig^2*xf_xf + 4/24*out.hsssv*sig^3*xf + 1/24*out.hssss*sig^4;
    
                %order 2
                xf_xf_p           = kron(xf_p,xf_p);

                %order 3
                xf_xs_p           = kron(xf_p,xs_p);
                xf_xf_xf_p        = kron(xf_p,xf_xf_p);

                %order 4
                xf_xf_xf_xf_p     = kron(xf_p,xf_xf_xf_p);
                xf_xrd_p          = kron(xf_p,xrd_p);
                xf_xf_xs_p        = kron(xf_xf_p,xs_p);
                xs_xs_p           = kron(xs_p,xs_p);
    
                Y_sim_1st(:,t)    = out.gv*xf_p;
                Y_sim_2nd(:,t)    = out.gv*xs_p + GGxxtil*xf_xf_p + 1/2*sig^2*out.gss; 
                Y_sim_3rd(:,t)    = out.gv*xrd_p + GGxxtil*2*xf_xs_p + GGxxxtil*xf_xf_xf_p + 3/6*out.gssv*sig^2*xf_p + 1/6*sig^3*out.gsss;
                Y_sim_4th(:,t)    = out.gv*xth_p + GGxxtil*(2*xf_xrd_p+xs_xs_p) + GGxxxtil*(2*xf_xf_xs_p) +...
                                    GGxxxxtil*xf_xf_xf_xf_p + 3/6*out.gssv*sig^2*xs_p + GGssxxtil*sig^2*xf_xf_p + 4/24*out.gsssv*sig^3*xf_p + 1/24*sig^4*out.gssss;
    
                X_sim_1st(:,t-1)  = xf_p;
                X_sim_2nd(:,t-1)  = xs_p;
                X_sim_3rd(:,t-1)  = xrd_p;    
                X_sim_4th(:,t-1)  = xth_p;    
           
                xf                = xf_p;
                xs                = xs_p;     
                xrd               = xrd_p;        
                xth               = xth_p;

                %oder 2
                xf_xf             = xf_xf_p;
    
                %order 3
                xf_xs             = xf_xs_p;
                xf_xf_xf          = xf_xf_xf_p;
    
                %order 4
                xf_xf_xf_xf       = xf_xf_xf_xf_p;
                xf_xrd            = xf_xrd_p;
                xf_xf_xs          = xf_xf_xs_p;
                xs_xs             = xs_xs_p;

            end

            Y_sim       = Y_sim_1st + Y_sim_2nd*(setup_.order>1) + Y_sim_3rd*(setup_.order>2) + Y_sim_4th*(setup_.order>3);
            X_sim       = X_sim_1st + X_sim_2nd*(setup_.order>1) + X_sim_3rd*(setup_.order>2) + X_sim_4th*(setup_.order>3);

            sim_total_t = vertcat(Y_sim,X_sim(1:end-M_.exo_nbr,:));

            names_dynplus    = vertcat(out.label_y,out.label_v(1:end-M_.exo_nbr,:));  

            for i_index=1:M_.endo_nbr
                name_temp                = setup_.names(i_index,:);
                sim_total(i_index,:)     = ys(i_index,:) + sim_total_t(strmatch(name_temp,names_dynplus,'exact'),:); %#ok<*MATCH3>
            end
        end
    end
    
    
        simulstore(1 + (s_index - 1)*endosize: endosize + (s_index - 1)*endosize,:)= sim_total; %save simulstore simulstore setup_.drop shocks_matrix names

    
    if setup_.numbersimul>1
       if (s_index == 1);                               disp(['s_index toolbox = ' num2str(s_index,'%.0f'),': Executed  0 percent of loop iteration']); end         
       if (s_index == round(.2*setup_.numbersimul));    disp(['s_index toolbox= '  num2str(s_index,'%.0f'),': Executed  20 percent of loop iteration']); end         
       if (s_index == round(.4*setup_.numbersimul));    disp(['s_index toolbox= '  num2str(s_index,'%.0f'),': Executed  40 percent of loop iteration']); end         
       if (s_index == round(.6*setup_.numbersimul));    disp(['s_index toolbox= '  num2str(s_index,'%.0f'),': Executed  60 percent of loop iteration']); end         
       if (s_index == round(.8*setup_.numbersimul));    disp(['s_index toolbox= '  num2str(s_index,'%.0f'),': Executed  80 percent of loop iteration']); end         
    end
end

output_toolbox_.simulstore = simulstore;    


if setup_.gen_irf==1
    irf_dev_matrix                = zeros(endosize,setup_.lengthsimul);
    ergodic_mean                  = zeros(endosize,1);
    for i_index=1:endosize
        x_noimpulse               = simulstore(i_index : endosize : endosize*setup_.numbersimul/2, :);
        x_impulse                 = simulstore(i_index + endosize*setup_.numbersimul/2 : endosize : endosize*setup_.numbersimul, :);
        sim_diff                  = x_impulse - x_noimpulse;
        
        if setup_.numbersimul>2
        sim_diff_mean             = mean(sim_diff);
        else
        sim_diff_mean             = sim_diff;
        end
        
        irf_dev_matrix(i_index,:) = sim_diff_mean;
        ergodic_mean(i_index,1)   = mean(mean(x_noimpulse(:,setup_.drop:setup_.lengthsimul)));
        if (max(max(abs(sim_diff(:,1:setup_.drop-1))))~=0); disp('mistake in dynare toolbox at setup_.gen_irf = 1: dynare_toolbox'); open dynare_toolbox; keyboard; end
    end
    
    yss                                = ergodic_mean;
    irf_dev_matrix(:,1:setup_.drop-1)  = [];
    irf_matrix                         = irf_dev_matrix;
end


if setup_.stoch_ss_irf==1
    irf_matrix              = simulstore;
    yss                     = irf_matrix(:,setup_.drop-1);
    irf_dev_matrix          = irf_matrix - repmat(yss,[1 setup_.lengthsimul]);
    
    if options_.periods==0
    irf_dev_matrix(:,1:setup_.drop-1)  = [];
    end
    
     irf_matrix_growth = irf_matrix - lagmatrix(irf_matrix',1)';
     
    
    if options_.periods<1 && max(max(abs(irf_matrix_growth(:,setup_.drop-1))))>1e-04
        disp('Non-stochastic ss has NOT been achieved');
        if setup_.nokeyboard==0
            keyboard
            open dynare_toolbox;
        end
    end
    
    if options_.periods<1 && max(max(abs(irf_matrix_growth(:,1:setup_.drop-1))))<1e-06 && setup_.order>1
        disp('No dynamics to reach the non-stochastic ss');
        %keyboard
        open dynare_toolbox;
    end
    
end

if (setup_.gen_irf==1 || setup_.stoch_ss_irf==1)
    output_toolbox_.irf_matrix     = irf_matrix; 
    output_toolbox_.irf_dev_matrix = irf_dev_matrix;
    output_toolbox_.yss            = yss;
    output_toolbox_.simulstore     = simulstore;
    output_toolbox_.shocks_matrix  = shocks_matrix;
    
    if setup_.gen_irf==1
        output_toolbox_.sim_diff       = sim_diff;

    end    
end

